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ABSTRACT 

J3 ■ We report evidence for extended gamma-ray emission from the Virgo, Fornax and Coma 

clusters based on a maximum-likelihood analysis of the 3-year Fermi-LAT data. For all 
three clusters, excess emission is observed within three degrees of the center, peaking 
at the GeV scale. This emission cannot be accounted for by known Fermi sources or by 
the galactic and extragalactic backgrounds. If interpreted as annihilation emission from 
super symmetric dark matter (DM) particles, the data prefer models with a particle 
■ mass in the range 20 — 60 GeV annihilating into the bb channel, or 2 — 10 GeV and 

> 1 TeV annihilating into fJ,~^ fj,^ final states. Our results are consistent with those 
obtained by Hooper and Linden from a recent analysis of Fermi-LAT data in the region 
of the Galactic Centre. An extended DM annihilation profile dominated by emission 
from substructures is preferred over a simple point source model. The significance of 
DM detection is 4.4(T in Virgo and lower in the other two clusters. We also consider the 
possibility that the excess emission arises from cosmic ray (CR) induced gamma-rays, 
and infer a CR level within a factor of three of that expected from analytical models. 
However, the significance of a CR component is lower than the significance of a DM 
component, and there is no need for such a CR component in the presence of a DM 
component in the preferred DM mass range. We also set flux and cross-section upper 
limits for DM annihilation into the bb and channels in all three clusters. 
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Introduction 



The existence of dark matter (DM) in the universe has so far only been deduced and constrained 
from its gravitational effect, due to the lack of electromagnetic interactions of the DM with itself 
or with baryonic matter. There are several elementary particle candidates for DM in various 
extension of the standard model of particle physics. Weakly interacting massive particles (WIMPs) 
are one of the most promising classes of dark matter candidates, with a self-interaction cross-section 
estimated at around 3 x 10~^^cm^s~^ under simplified assumptions. These particles arise naturally 
in theories that seek to extend the standard model and, at the same time, produce the correct relic 
density of dark matter in the early universe. Within the framework of the minimal supersymmetric 
standard model (MSSM), the lightest neutralino emerges as the prototype of a WIMP that is stable 
over cosmological timescales and can annihilate into standard model particles. These particles are 
generally known as cold dark matter. 

Much effort has been devoted to the search for dark matter, either directly or indirectly. Direct 
detection involves identifying the rare events of DM scattering off ordinary matter or searching for 
new particles near the weak scale at the LHC. Indirect detection involves looking for the annihilation 
or decay products of dark matter in cosmic rays and gamma rays. In particular, pair annihilation 
produces gamma-ray photons at a rate proportional to the square of the dark matter density, 
which then propagate, almost without absorption, to the o bserver. In this cas e, the Galactic 
Center should be the brightest gamma-ray source on the sky (jSpringel et al.l 120081 . and references 
therein). There have b een recent claims of extended emission (besides a central point source 



2011 



in the Galactic Center (IHooper &: LindenI l201ll : iHooper &: Goodenoughl l201ll : IChernvakova et al 



Boyarsky et al.ll2010l ). which is consistent with either DM annihilation from a halo with a 
slightly cuspier densi ty profile (an inner slope o f ~ —1.25 to —1.4) than the typical NFW cold dark 



matter galactic halo (jNavarro et al 



1996 



19971 ) , or with proton collisions accelerated by the central 



galactic black hole. The latter process is very poorly understood so a definitive conclusion is not 
yet possible. 

Targetting the entire sky rather than the Galactic center in searching for annihilation radiation 
may seem a good strategy since it takes advantage of the large-scale distribution of dark matter while 
avoiding some of the uncertainties arising from the astrophysical modelling of galactic gamma-ray 
sources. However, the fact that we are located near the center of the Galactic halo and that most 
of the annihilation emission outside the Ga lactic Center is produced by dark matter substructures 



(jSpringel et al, 



2008 



Diemand et al.l 120071 ) results in a gamma-ray map from annihilation that is 



almost uniform on large scales. This makes detection within the Milky Way halo a difficult task, 
exacerbated by the additional unc ertainty in having to rnodel the extragalactic ba ckground, which 



is more important on large scales (jZaharijas et al.ll2010l : [Baxter fc Dodelson 



2011 



Dwarf galaxies are the most DM-dominated objects known, are relatively free from astrophysi- 
cal contamination and appear compact on the sky. They are therefore p romising targets for DM an- 

nihilation detection. Recent joint analyses of eight to ten dwarf galaxies ()Geringer-Sameth &: Koushiappas 



- 3 - 



201 ll : iThe Fermi-LAT Collaboration: M. Ackermann et alj 120111 ) resulted in a null detection but 
have began to rule out the canonical annihilation cross-section of 3 x 10~^^cm^s~^ for DM masses 
below ~ 20 GeV. 

Galaxy clusters are the most massive virialized DM structures in the universe and are also 
good targets for indirect DM searches. The presence of a large population of DM substructures (or 
subhalos) predicted by numerical simulations further boosts the detectability of DM in clusters. 
Although the total mass within subhalos amounts to only 10 to 20 percent of the total halo mass, 
the density enhancement within subhalos can boost the total cluster annihilation luminosity by 
a factor as high as 1000 when extrapolated down to a subha. l o ma s s limit of one Eart h mass as 
expected, for example, for a ~ 100 GeV WIMP (|Gao et al.1 boill : IPinzke et aPboill ). As the 



distribution of subhalos is much less concentrated than that of the smooth main halo, the total 
annihilation emission from clusters is predicted to be extended. Thus, attempts to detect DM 
annihilation assuming a point source or NFW-squared type profile could miss most of the signal. 
In fact, j ust such a search using the 11-month Fermi-LAT data has yielded a null detection in six 
clusters (j Ackermann et al.ll2010l ) 



In this work, we use the 3-year Fermi-LAT data to search for extended emission from clus- 
ters. We consider possible contributions from cosmic ray (CR) induced g amma-ray emission 
which can be as high as, or higher than, tha t from cluster DM annihilation (IJe tema et a 



Pinzke &: Pfrommer 



2010 



Pinzke et al. 



2011 



We adopt the model proposed by 



Gao et al 



2009': 
(1201 ih 



for th e cluster annihilation profile and the semi-analytic method developed by lPinzke &: Pfrommer 
(|2010l ) to model CR induced gamma-ray emission and provide constraints on both the DM and 
CR components. We focus on three ga l axy c lusters: Coma, which is predicted to have the highest 
signal-to-noise according to iGao et al.l (120111). and Fornax and Virgo which are predicted to have 



the lowest astrophysical contamination (jPinzke et al. 



2011 



During the final stages of preparation of this work, a paper (jHuang et al.ll201ll ) was posted 
on arxiv presenting a null detection of DM annihilation emission from a combined analysis of 
eight galaxy clusters. Our work differs from this interesting paper in several respects: firstly, we 



assume a DM annihilation profile based on high resolution cosmological simulations (jGao et al 



2011 



); secondly, we assess the impact of cosmic rays in the detection of dark matter; and finally, we 
include the Virgo cluster in our sample which turns out to be the b est candidate . The constraints 
we set on the annihilation cross-section are consistent with those in lHuang et al.l (|201lh . 



(2011 



The cosmological parameters used in this work are the same as those assumed by lGao et al 



0.25,J^A = 0.75,/i = 0.73. 



2. Modeling gamma-ray emission in clusters 



We model the observed gamma-ray emission in clusters with several components as shown in 
Figured) the galactic foreground (GAL), the extragalactic background (EG), emission from known 
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Fig. 1. — Decomposition of the Fermi-LAT image in the region of the Virgo cluster into model 
components. The observed photon count image from lOOMeV to 100 GeV is shown on the left. 
The right panels show the integrated image over the same energy range for the various model 
components: galactic diffuse emission, extragalactic diffuse emission, 2FGL point sources, cosmic- 
ray photons and DM annihilation emission, as labeled. The green circle in the "Point Sources" 
panel marks the virial radius of the cluster. The "Fluctuation" panel shows the residual image for 
our best-fit DM model. The images have been enhanced individually in color space for contrast. 



point sources, DM annihilation and CR-induced emission. The GAL and EG diffuse emission are 
given by the most recent templates, gal_2yearp7v6_v0 . f its and iso_p7v6source.txt, which can 
be obtained from the Fermi-L AT data server, while the point s o urces are taken from the LAT 2- year 
point source catalogue (2FGL jThe Fermi-LAT CollaborationI (j201ll )). We now describe in detail 
our models for DM annihilation and CR emission. 



2.1. Cluster annihilation emission 

The gamma-ray intensity along the line-of-sight due to DM annihilation is given by: 

dNf 

where is the DM particle mass, the particle model dependent term giving the differential 
number of photons produced from each annihilation event as a function of energy, ii^, in a particular 
annihilation channel, /, and afV is the cross-section (or annihilation rate) for that channel, which is 
predicted to be constant for cold dark matter. The line-of-sight integration of the density squared 
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is often expressed in terms of a dimensionless factor, 

If the source size is much smaller than the instrumental beam size, a point source approximation 
is applicable. In this case, the integration of J over a large enough solid angle is used to determine 
the total flux for the point source, Jint = /^^ JdVL. 

The cluster annihilation emi ssion is modeled w ith both a point source approximation and with 



the extended profile suggested by lGao et al.l (j201ll ). We will refer to models using these two profiles 



respectively as PT and EXT. If the cluster follows a smooth NFW profile, then its integrated J 
factor which determines the total annihilation flux can be found as 

47r 2 3 1 1 / 1 \2 

Jnfw = -jPsTs^^ X ^(o,3GeV/cm3) " 

Here Da is the angular diameter distance of the cluster and ps and are the characteristic density 
and radius for the NFW profile. They are related to halo concentration and virial radius through 

the relations, = — -, ^ and Vg = rono/c, with pr the critical density of the uni- 

3 log(l + c) - c/(l + c) ' ' ^ ^ 

verse, r2oo the cluster virial radius within which the average density is 200pc and the concentration 
parameter, c, given by the following mass-concentration relation: 

c = 5 74( ^ (4) 



(jPuffv et al.ll2008l ). Extrapolating to a cutoff mass of 10 ^Mq, the existence of subhalos will boost 



this flux by a factor 

b{M2m) = Jsuh/JNFW = 1.6 X lO-3(M2oo/M0)O-39 (5) 



Gao et al.l (j201ll ). Using the results of the simulations by these authors, the surface brightness 



profile of subhalo emission can be fitted within r2oo by the following formula: 

, , 166(M2oo) Jatfvk 1 / / \ tci\ 

.ln(17) r|oo + 16r2 

Below we fit the subhalo emission surface brightness beyond the virial radius and extrapolate to 
several times the virial radius using an exponential decay, 

Jsubir) = J.n6(r20o)e-2-^^^(''/^-»-^) (r > r2oo). (7) 

The total annihilation profile is the sum of the contributions from a smooth NFW profile and the 
subhalo emission. This is completely dominated by subhalo emission except in the very center 
of the cluster. We show the total annihilation profile and its decomposition into main halo and 
subhalo contributions in the left panel of Figure [3l taking Virgo as an example. This profile is 
further inflated after convolution with the LAT point spread function. 
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Following I Ackermann et alj (|20ld ). we consider two representative annihilation channels, namely 
into b — b and — m~ final st ates. The annihilation spectrum is calculated using the DarkSUSY 
package ( Gondolo Sz Silk 19991 )1^ which tabulates simulation results from PYTHIA0 We also in- 
clude the contribution from inverse Compton (IC) scattered photons by energetic elec tron-positron 
pairs produced during the annihilation process, following the procedure described in lPinzke et al 



(|201ll ) . In general three external energy sources are involved in the dissipation and scattering of the 



injected electrons from annihilation: the Cosmic Microwave Background (CMB), in frared to UV 



light f rom stars and dust, and the interstellar magnetic field. However, as shown in iPinzke et al 



(|201ll ). the latter two components are expected to be important only in the inner region of clusters 
(< 0.03i?200c)) corresponding to less than 0.1 degrees for our three clusters. Including them would 
introduce a position-dependent component to the annihilation spectrum, so for simplicity we only 
consider the contribution of CMB photons in the IC calculation. For the bb channel, IC photons 
only contribute significantly to the low energy spectrum for relatively high neutralino mass, while 
for the ^"'"/U" channel, which has plenty of energetic electrons, the IC emission can completely 
dominate the annihilation emission over the full energy range of interest for the highest neutralino 
masses considered. 



We note that the electroweak corrections recently proposed bv lCiafaloni et al.l (|201lh (see also 
Cirelli et al.l (|201ll )) can bring visible differences to the channel spectrum at high WIMP 

masses before IC scattering. However, since IC photons dominate at the high mass end and the 
electroweak correction only significantly changes the positron yields at low energy, thus having 
little effect on the IC spectrum, the electroweak correction to the total spectrum is still negligible. 
The total photon yields are shown in Figure El The almost flat spectrum with a cutoff around 
the energy corresponding to the WIMP mass comes from prompt annihilation emission including 
continuum secondary photons and final state radiation from charged final state particles. The low 
energy rise originates from IC scattered CMB photons. 



2.2. Cosmic-ray induced gamma-ray emission within clusters 



The cosmic ray induced gamma-ray emission is calculated following a semi-analytic prescrip- 
tion, derived from high r esolution numerical simulat ions of galaxy clusters, that models cosmic ray 
physics self consistently (IPinzke &: Pfrommeii 12010 ). The gamma-ray photon production rate (or 
source function) from pion decay is found to be separable into a spatial and a spectral part: 



qcR{r,E) = 



dtdVdE 



A[r)s{E), 



(8) 



""^http: / /www.darksusy.org. 

^http: //home. thep.lu.se/ torbjorn/Pythia.html 
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Fig. 2. — Photon yields for bb (red lines) and fJL^ (green lines) channels. Plotted are the total 
photon yields including continuum secondary photons, final state radiation from charged final state 
particles, as well as inverse Compton scattering of CMB photons by electron/positron pairs, for the 
mass range 10— 1000 GeV for the bb channel and IGeV — lOTeV for the channel. The masses 

are sampled uniformly in a log scale. Note that each spectrum cuts off at an energy corresponding 
to the particle mass. For comparison, the black line shows the photon spectrum from cosmic ray 
induced photons with arbitrary normalization. 



- 8 - 



where the spatial part, A{r), is proportional to the square of the gas density profile multiplied by 
a slowly varying radial function parametrized by cluster mass. The spectral part, s{E), is almost 
independent of cluster mass and has a power-law form, dN^ / d\n.{E^) cx E~^'^, for the energy range 
1 ~ 100 GeV but flattens at low energies, as shown in Figure [2j We summarize the detailed form 
of A{r) and s{E) plus the gas density profile for the three clusters derived from X-ray observations 
in the Appendix. 

The differential gamma-ray flux from this source function is simply its integral along the line- 
of-sight: 



This prescription is derived from the average emission profile for a sample of simulated clusters 
for a realistic choice of parameter values (e.g., for the maximum shock acceleration efficiency, 
Cp,max)- In addition to the uncertainties in the model parameters there is also uncertainty in the 
observationally derived halo mass and gas density profile. In this work, we simply assume that the 
shape of qcR{r,E) is given by the model described above and account for the uncertainty in the 
model parameters as well as sample variance with an additional normalization parameter, acR, so 
that 



We take acR = 1 as our fiducial CR model and also consider the case when acR is fitted from 
the actual gamma-ray data as an optimal model. In the right panel of Figure [3] we compare the 
CR profile for the fiducial model to the expected DM annihilation profile within our three clusters, 
assuming a fiducial DM particle model with particle mass M ~ lOOGeV annihilating through the 
bb channel with cross-section, av = 3 x 10~^^cm^s~^. In general the CR emission is more centrally 
concentrated than the annihilation profile since the CR traces the gas profile. It can be readily 
seen that Fornax has a particularly low CR level while Coma is CR dominated. Coma has steeper 
profiles due to its larger distance and hence smaller angular size. 




(9) 




(10) 
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Fig. 3. — Cluster photon profiles. Left: theoretical and PSF-convolved J profile for Virgo. The 
total annihilation profile is shown as a black solid line and is decomposed into the smooth main 
halo part (red dashed line) and the subhalo part (blue dashed line). The green solid line shows the 
annihilation profile after PSF convolution, plotted down to an inner radius comparable to the pixel 
size of 0.1 deg. Right: PSF-convolved photon profiles from annihilation (solid) and cosmic rays 
(dashed) for three clusters (indicated by different colors). Solid lines show the expected photon 
count profile for the extended DM annihilation model. Dashed lines show the expected cosmic-ray 
induced photon counts for the fiducial CR model. For comparison, we also plot the PSF-convolved 
profile for a central point source model (black solid line) with arbitrary normalization. In both 
panels, a dark matter model with particle mass, M ~ lOOGeV, and annihilation cross-section, 
av = 3 X lO^^^cm^s^^, through the bb channel is assumed. 
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3. Data analysis 

3.1. Data preparation 

We analyze the first 3 years of data (04/08/2008 to 16/08/2011) from the Fermi-LAT, 1 
trimmed with the cuts listed below, to select high quality photon events. This typically results 
in ~ 10^ photons within a radius of 10 degrees around each cluster, while the expected number 
of annihilation photons is of the order of lO'^ according to Figure [3j The most recent instrument 
response function, P7SOURCE_V6, is adopted for the analysis, in accordance with our event class 
selection. 



Minimum Energy 100 MeV 

Maximum Energy 100 GeV 

Maximum zenith angl^ 100 degrees 

Event ClasJl 2 (P7S0URCE) 

DATA-QUAlfl 1 

LAT CONFiqll 1 

ABS (ROCK ANGLeJ < 52 

ROI-based zenith angle cut yes 



We list the basic properties of the three clusters in Table [TJ 



■^http: / /fermi. gsfc.nasa.gov /cgi-bin /ssc /LAT /L ATDataQuery.cgi 

■'zenith angle (degrees): angle between the reconstructed event direction and the zenith line (originates at 
the center of the Earth and passes through the center of mass of the spacecraft, pointing outward). The Earth's limb 
lies at a zenith angle of 113 degrees. 

^EVENT CLASS: flag indicating the probability of the event being a photon and the quality of the event recon- 
struction. 

®DATA-QUAL: flag indicating the quality of the LAT data, where 1 — OK, 2 — waiting review, 3 — good with 
bad parts, — bad 

^LAT-CONFIG: flag for the conflguration of the lat (1 = nominal science configuration, = not recommended 
for analysis) 

*ROCK ANGLE: angle of the spacecraft z-axis from the zenith (positive values indicate a rock toward the north). 



- 11 - 



Table 1: Basic Properties of Target Clusters 





Coma 


Fornax 


Virgo (M87) 


RA (deg) 


194.9468 


54.6686 


187.6958 


DEC (deg) 


27.9388 


-35.3103 


12.3369 


Da (Mpc)^ 


95.8 


17.5 


22.3 


M200 (Mo)^ 


1.3el5 


2.4el4 


6.9el4 


^^200 (deg)^ 


1.3 


4.1 


4.6 


Jnfw^ 


5.9e-5 


4.1e-4 


6.3e-4 


Boost factor'^ 


1.3e3 


6.5e2 


9.8e2 



"Angular diameter distance from the NASA extragalactic database. 

''Cluster halo mass defined as mass within the radius, r2oo, within which the a verage density equals 200 times the 
critical density of the universe. Values of masses taken from lPinzke et al.l (|201ll ). 

'^Integrated coefficient, jTint = /^^ JdQ, over the solid angle spanned by the cluster virial radius, assuming a smooth 
NFW density profile. 

''Boost factor relative to the total annihilation luminosity within the virial radius due to substructures, extrapolated 
to a subhalo mass limit of 10"^Mb 



3.2. Maximum-likelihood fitting 

We use the pyLikelihood tool shipped with the Fermi Science Tool s software package ( version 



v9r23pl-fssc-20111006) to perform a maximum likelihood (ML) analysis (jMattox et al.lll996l ). After 
applying appropriate data cuts, as described in section [3TT| we bin the data into 0.1 degree-wide 
pixels and 30 logarithmic energy bins within a radius of 10 degrees around each cluster. This large 
radius is chosen to account for the large LAT PSF size (4 ~ 10 degrees at 100 Mev|§). 

A model is constructed to fit the data including all known foreground and background emis- 
sion, as well as DM and CR components. We include all the point sources from 2FGL within a 
radius of 15 degrees from the cluster center in the model, plus the most recent galactic (GAL) 
and extragalactic (EG) diffuse emission given by the template files gal_2yearp7v6_v0 . f its and 
iso_p7v6source.txt. The normalization of the GAL and EG diffuse components are allowed to 
vary during the fitting. Within the cluster virial radius there are two 2FGL point sources in Fo rnax, 



and three in Virgo, one of which is associated with the central AGN, M87 (jAbdo et al.ll2009l ). We 
allow the normalization and power-law spectral index of these five point sources to vary freely. 
Parameters for other point sources are fixed as in the 2FGL catalog. 

A surface brightness template given by the dimensionless factor J in Equation [2] is generated 
for each cluster out to a 15 degree radius by summing up both the contribution from a smooth NFW 
profile and the boost from subhalos. This J-map is used to fit for extended cluster annihilation 



'The LAT PSF size scales roughly as E~°-^, so at 1 GeV it is ~ Ideg 
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emission. For the point source model, the integrated factor Jnfw in Equation[3]is used to derive an 
annihilation cro.;secfon from the attedtotal flux§ Similarly, a CR photon template is generated 
for each cluster out to three times the cluster virial radius, where the surface brightness has dropped 
to below 10~^ of the central value and beyond which the model is not reliable. Images for various 
model components are shown in Figure [1] taking Virgo as an example. 

The photon counts within each pixel are treated assuming Poisson statistics for each energy 
bin to calculate the likelihood. The best-fit parameters are obtained when the likelihood for the 
entire data set is maximized. The significance of a given component of interest (e.g. DM or CR) 
from the ML fitting is quantified by the likelihood ratio statistic, 



r5 = -21n(Lo/L), (11) 

where L is the maximum likelihood for the full model, and Lq is the maximum likelihood for the null 
hypothesis, i.e, the model without the component of interest. This test statistic (TS) approximately 
follows a distribution, with one degree of freedom for our case where the normalization is the 
only free parameter. The probability that a given value of TS arises purely from fluctuations of 
the null hypothesis is: 

P= / i;xl{m= / ^^dx. (12) 



TS 2 Jy/TS v27r 



The factor — comes from the constraint that the normalization parameter is non-negative. The 
significance of detection can thus be quoted as VTSa (one sided Gaussian confidence). 



4. Results 

4.1. Constraints on CR emission 

With all the model components defined above, we first proceed with ML fitting for a model 
with no DM annihilation but with cosmic rays, the "CR-only" model hereafter. Note that the 
GAL and EG backgrounds, as well as the nearby point sources are always implicitly included in 
the analysis, as described in section 13.21 The results for the CR-only model fitting are listed in 
Table [2j The fitted CR levels all agree within a factor of three with the theoretical predictions. 
While Fornax is the most poorly constrained due to its intrinsically low CR level, the derived upper 
limit for Coma already rules out the fiducial value at 95% confidence. 



^ This only strictly applies if the cluster size is smaller than the PSF size so that the entire cluster within the virial 
radius can be approximated as a point source. Otherwise, the integration over solid angle should only extend to the 
angular resolution of the telescope. However, since the surface brightness for a smooth NFW profile drops rapidly 
with radius, and because of the energy dependent PSF size, we choose to use a generic estimation of Jint as given by 
an integral over the cluster virial radius. 
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Table 2: Fitting to the CR-only Model 







0!CR,UL 


Fcr,ul'' (pli • g-i) 


TS 


Coma 


0.3 ±0.1 


0.5 


2.74e-09 


5.2 


Fornax 


2.0 ±2.6 


6.4 


2.4e-09 


0.6 


Virgo 


1.5 ±0.5 


2.3 


2.1e-08 


8.2 



"Best fit normalization {acR,fit = 1 is the theoretical prediction.) 
''95% upper limit (UL) on the normalization 

°95% upper limit on the CR induced gamma-ray flux from 100 MeV to 100 GeV 



4.2. Constraints on DM annihilation 

Given the low significance of the CR detection in the CR-only model, it is not safe simply to 
adopt the best fit acR values for further extraction of the DM signal. Instead, we consider the 
following four families of cosmic ray models in the presence of a DM component: 

Fiducial-CR model. The CR level is fixed to the theoretical expectation, acR = 1- Since this is 
larger than the upper limits derived for Coma, we exclude Coma from further discussion for 
this family. 

Optimal-CR model. The CR level is taken as the best-fit value listed in Table [2j 
FVee-CR model. The normalization of the CR level is treated as a free parameter. 
No-CR model. No CR emission is considered, only DM. 

For each family, both point source (PT) and extended (EXT) profiles are considered for the DM 
component. 

4.3. The 66 channel 

The significance for the 66 channel of annihilating DM is shown in Figure [H It can be seen 
that in almost all cases, the EXT profile has higher significance than the PT profile, over the entire 
DM mass range considered. There is no significant evidence for a PT model component within 
Fornax. For the four families of models, the No-CR model has the highest DM significance. In 
particular, a peak TS value of 19.8, corresponding to a 4. 4a confidence level, is obtained for Virgo 
with the EXT profile when no CR component is considered. This decreases to 3.4cr when a CR 



- 14 - 



component with free normalization is included and is even lower when the CR level is fixed either 
to the fiducial or the optimal levels. There is a general trend for the significance t o peak in the mass 
range 20 — 60 GeV for all three clusters. This is consistent with the conclusion of lHooper Linden 
(|201ll ) , who claim that a DM model with particle mass in the range 25 — 45 GeV annihilating into 
b b final states can explain the excess extended emission observed from the direction of the Galactic 
center. 




H10 



10 

M^/GeV 



M^/GeV 



Fig. 4. — TS values for a DM component annihilating through the bb channel. Different colors rep- 
resent different clusters. Solid lines correspond to the extended model while dashed lines correspond 
to the point source model. The four panels are for the various CR models as labeled. 



To facilitate comparison between all different model families and assess the significance of a 
CR component, we show the maximum likelihood values for each model of the Virgo cluster in 
Figure [5j As seen in the TS curves, the likelihood for the extended model is always higher than 
that for the point source model. For the most likely mass range, 20 ~ 60 GeV, the No-CR and 
Free-CR models with extended DM emission share the highest likelihood, with the No-CR model 
being superior by having one fewer parameter. Actually there is effectively no contribution from 
CR when CR and DM are fitted simultaneously for this particle mass range, and the fiducial and 
optimal CR levels are above the 95% CR upper limit predicted from the Free-CR model in the 
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Fig. 5. — Likelihoods for various model fits to the Virgo cluster data. Plotted is the logarithm of 
the maximum likelihood values relative to that for the no-CR model with no DM component, as a 
function of the DM particle mass. The colors correspond to the four families of CR models. For 
each model family, the solid line shows the likelihood when the EXT profile is assumed, while the 
dashed line shows the likelihood for the PT approximation. The dotted lines show the likelihood 
for models with no annihilation emission, in which case the Free-CR and Optimal-CR models are 
equivalent. 
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presence of extended DM annihilation. This preference for DM over CR also obtains for Fornax, 
although the optimal and fiducial CR cases are still allowed in the presence of an extended DM 
component. The story is very different in Coma, where the DM and CR components are almost 
degenerate in the sense that either component alone has a similar level of likelihood. 

In Figure [6] we show the 95% confidence upper limits on the DM annihilation flux and compare 
them to the CR levels. For each cluster, the colored stripes are defined by the minimum and 
maximum upper limit predictions among the four families of models. The optimal CR levels in 
the three clusters are all comparable to the fitted DM flux, and the DM flux upper limits for the 
four different CR models vary only within a factor of two, with the No-CR and Optimal-CR cases 
predicting the highest and lowest upper limits. The left and right panels show the results for the 
PT and EXT models respectively; the PT models always have lower flux upper limits than the 
extended models. 

The flux upper limits are translated into cross-section upper limits in Figure [71 using Equa- 
tion [TJ These are also shown as colored regions reflecting the variance in the different treatments of 
CR. Although the predicted flux upper limits decrease slowly with DM particle mass and remain 
within the same order of magnitude for the mass range considered, the resulting cross-section upper 
limits increase by a factor of 100 from low to high particle mass. This is because low mass particles 
correspond to higher DM number densities (the p^/M^ factor in Equation [T]) for a given mass 
density, so to obtain the same flux level, the required cross-section must be smaller for low mass 
particles. With the help of a boost factor of order 10^, a much lower cross-section is needed (by a 
factor of at least 100) for extended annihilation models to achieve a slightly higher flux upper limit 
than point source models. 



Our cross-section limits are much lower than the 11-month Fermi-LAT constraints in Ackermann et al 



(|2010l ). where the tightest constraint came from Fornax for a much lower assumed boost factor of 
~ 10. Our lim its are also tighter than that frq n i a jo int analysis of the dwarf satellites of the 



Milky Way by iGeringer-Sameth fc Koushiappaa (120111 ). They drop below the fiducial thermal 



cross-section of 3 x 10~^^cm^s~^ for < lOOGeV. Of the three clusters, Virgo has the highest sig- 
nificance and flux upper limits and places the tightest constraints on the annihilation cross-section. 
The best-fit cross-section for the extended model is closer to the thermal cross-section than that for 
the point source model. For example, for M-^ = 28 GeV and where the TS value is higher than 2, 
the best-fit cross-section is (3 — 6) x lO^^'^cm'^s^^ for the extended model and (2 — 4) x 10~^^cm^s^^ 
for the point source model. 



4.4. The channel 

The significance, flux upper limits and cross-section upper limits for DM annihilating through 
the channel are plotted in Figures [8] to [TOl The predicted flux upper limits for Coma and 

Virgo are still comparable to the CR level, with Fornax having much lower CR emission. A DM 
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Fig. 6. — DM annihilation flux upper limit for the bb channel. The stripes are defined by the 
minimum and maximum upper limits given by the four CR model families, with different colors 
corresponding to different clusters, as indicated in the legend. Left and right panels are the results 
for the EXT and PT profiles respectively. For each cluster, a solid line of the corresponding color 
shows the optimal CR flux. 
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Fig. 7. — Upper limit for the DM annihilation cross-section in the bb channel. The different 
colors represent the three clusters, with the stripes spanning the range between the minimum and 
maximum upper limits given by the four different ways of treating the CR component. The three 
highest stripes show the PT model constraints and the three lowest the EXT model constraints. 
We also plot with dashed lines constrai ntj^ from a joint analysis of the Milky Way dwarf galaxies 
( Geringer-Sameth Sz Koushiappag 201ll. black dashed l i ne) a nd previous constraints from the 11- 
month Fermi-LAT data for Fornax (lAckermann et al.l l201Cll . green dashed line) assuming these 
authors' optimistic value for the boost factor from subhalos, which gives the tightest constraint. 
The black solid line indicates the canonical thermal cross-section of 3 x 10~^^cm^s~^. 



'^The lines for "Dwarf and "Fermi- lyr" are only reproduced schematically, by reading out 
several data points from the original plots in the references. 
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particle mass of 2 — 10 GeV is preferred, which is also cons i stent with the values inferred from 
analysis of the Galactic center emission bv lHooper &: LindenI (j201ll ). In addition, a second region 
of high significance is obtained for > 1 TeV. The inferred cross-section falls below the canonical 
value for DM particle masses less than 10 GeV. Note the discontinuity in the upper limit predictions 
around 100 GeV which reflects the transition from the prompt annihilation dominated regime to 
the IC emission dominated regime in the photon spectrum. 
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Fig. 8. — TS values for DM annihilating through the /U^/i channel. Line styles are as in Figured! 



4.5. Examination of the excess emission 

To pin down the source of the significance of our fits, in Figure [TT] we break down the TS value 
for our preferred particle model into contributions from four radial ranges. Clearly, most of the 
significance for the DM model comes from the region within 3 degrees around each cluster. The 
spectrum within this region (Figure [T2|) shows excess emission around 1 GeV for all three clusters 
relative to a fit with no DM or CR components (the null model). This excess is consistent with an 
almost scale-free spectrum (dN/dlnE ~ const), picking out DM particles in the range 20 — 60 GeV 
for bb and 2 — 10 GeV for final states, for which the cutoff at the high energy end and the IC 
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Fig. 9. — DM annihilation flux upper limits in the ^'^ n channel. Line styles are as in Figure El 



boost at the low energy end are not significant around GeV energy scales. At much higher masses 
(above 1 TeV) in the channel, where IC emission completely dominates the annihilation 

spectrum and begins to harden again around GeV scales, the TS for the annihilation component 
picks up again, as seen in Figure El 

The radial profile of the photon counts also reveals signs of excess emission in the central 
regions of the clusters when fitted with the null model, as shown in Figure [131 The central cluster 
emission is apparent in Coma and Fornax but in Virgo it can be observed only weakly in a ring of 
[1,3] degrees, due to the existence of a central point source, M87. We show the fitted photon counts 
spectrum of M87 in the null, DM-only and CR-only models explicitly in Figure [HI The total fiux 
for M87 is comparable to the CR and DM flux and varies by 20 to 30 percent in different models 
without much change in the spectral shape. The low energy dip for the three components comes 
from the decrease in the instrument's effective area at low energy. The AGN spectrum differs from 
the CR spectrum primarily at low energy and is distinct from the DM spectrum over a wide energy 
range. This, together with their different radial profiles, helps to disentangle the two components 
easily during the fitting. In Fornax, however, the central excess is only partly fitted by a DM 
component. We checked that the remaining part is mostly associated with a low energy spectral 
excess below 400MeV, which is not fitted by the current DM spectrum. 

A possible reason for a central excess we have found in all three clusters is a systematic error 
in the fitting. For instance, if the fixed point sources outside the cluster virial radius have too 
much flux attributed to them, then the diffuse EG and GAL components would be biased low and 
would not contribute enough model flux to the central region. This would lead to an excess in the 
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Fig. 10. — Upper limits for the DM annihilation cross-section in the u' ^ u~ channel. Line styles are 
as in Figure [71 The green dashed line is the 11-month Fermi result (lAckermann et al.ll201Cl ) for 
Fornax. 
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center. However, we have verified that our fits are insensitive to whether all the point sources are 
kept fixed or allowed to vary freelvl^ In addition, we have carried out the same analysis on three 
regions randomly selected from the 3-year LAT sky where no known prominent clusters exist near 
the center. No central excess is observed in any of these three regions, suggesting that the excess 
emission is related to the presence of the galaxy clusters. 

To see how common the extended emission is in galaxy clusters, we have searched in seven 
other clusters adopting our preferred DM particle model with ~ 28 GeV in the bb channel. 
These seven clusters ~ M49, A1367, AWM 7, Centaurus Hydr a, NGC4636 and NGC5813 - are 



predicted to have high annihilation flux by iPinzke et al.l (|201lh . They are also selected to lie at 



least 15 degrees away from the galactic plane. A DM significance comparable to that in Fornax is 
found in M49 and Al36£l with a TS value around 4, while Centaurus, NGC4636 and NGC5813 
show no sign of DM annihilation emission at all. We point out that these additional clusters are all 
very compact on the sky (except M49) with a virial radius close to or smaller than 2 degrees, making 
it difficult to exploit the information contained in the radial profile. In addition, the existence of 
additional sources inside the virial radius (one point source in Hydra and NGC4636, two point 
sources and one extended source (CenA Lobe) inside Centaurus) complicates the analysis further. 



^^When the M87 parameters are fixed to the 2FGL values, the significance of the DM detection in Virgo decreases 
to less than 3cr. However, we find that the 2FGL parameter values tend to overfit the current data slightly. When 
the M87 parameters are refitted within the null model, the normalization decreases by 15 percent and if this value is 
then fixed when fitting for an additional DM component, the significance for the DM and CR components are almost 
the same as when the fits are performed in our standard way. 

^■^The central excess in A1367 is more consistent with a point source profile, with a significance of 3.5a. 
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Fig. 11. — TS profile for our three clusters in the case of a DM model with ~ 30 GeV, 
annihilation through the hh channel and an extended emission profile. The TS value for each 
cluster is broken down into contributions from each of four radial ranges, [0,3], [3, 5], [5, 7] and [7,10] 
degrees, by calculating the likelihood ratio inside each region for the global best fit parameters. 
The colors correspond to the three clusters as labelled. The significance is dominated by regions 
within 3 degrees around each cluster center. 
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Fig. 12. — Observed and fitted energy spectrum in the inner 3 degrees of eacli cluster. The top 
panels show the counts (circles with error-bars) and the null model (the model without CR or DM; 
solid line) fit in each energy bin, multiplied by the energy of the bin to reduce the dynamic range. 
The lower panels show the residual for the no-CR model (black dashed line), the best-fit CR model 
without DM (green dashed) and the best-fit DM model without CR (red dashed) for the case of 
~ 30 GeV annihilating into the bb channel, normalized by the estimated Poisson error in each 
bin. For comparison we also show the contribution from the CR (green solid line) and DM (red 
solid line) components in the corresponding models. For Virgo the black solid line in the lower 
panel shows the contribution from the central AGN in the null model. The best-fit parameters 
are taken from the global best fit, i.e, from fitting the entire 10 degree region. Note that we have 
rebinned the data into 10 energy bins to produce this plot. 
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Fig. 13. — Observed and fitted radial photon profile around each cluster. The top panels show the 
counts (circles with error-bars) and the fits for three models: the null model (i.e. without CR or 
DM; black solid line), the DM-only model {My. ~ 30 GeV annihilating into the hh channel, red 
solid line) and the CR-only model (green solid line). The vertical dashed lines mark the cluster 
virial radii. The lower panels show the residuals for the three models, normalized by the estimated 
Poisson error in each bin. The best-fit parameters are taken from the global best fit, i.e, from 
fitting the entire 10 degree region. 
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Fig. 14. — Contribution of the central AGN in Virgo (M87) in different models. Circles with error- 
bars are the observed photon counts in each energy bin. The black solid line shows the best-fit 
model without a DM or CR component (the null model). The red solid line gives the counts for 
the extended DM component in the DM-only model and the green solid line gives the counts for 
the CR component in the CR-only model. The three dashed lines give the fitted counts for M87 
in the null (black), DM-only (red) and CR-only (green) models respectively. Only counts within a 
radius of three degrees are shown. 
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5. Discussion and conclusions 

We have performed maximum likelihood fits to the 3-year Fermi-LAT data for three galaxy 
clusters: Coma, Fornax and Virgo. We fit models which, in addition to point sources and galactic 
and extragalactic backgrounds, include emission due to dark matter (DM) annihilation and cosmic 
rays (CR). For the former, we assume both a point source and the theoretically predicted extended 
distribution of gamma rays in two generic annihilation channels. When searching for a dark matter 
signal, we experiment with different treatments of the CR component. Our main results are as 
follows. 



Assuming no DM annihilation radiation, the gamma ray data from Coma and Virgo already 
set significant constraints on the CR level at more than 2a. For Virgo, the data are co n sisten t 
with the prediction s from the analytic CR model proposed bv iPinzke fc: Pfrommerl ([20101); 
Pinzke et al.l (1201 ll ) while, for Coma, the data give an upper limit that is a factor of two lower 
than the analytic prediction, indicating either an uncertainty in model parameters such as 
halo mass and gas density, or a peculiarity of the CR emission in Coma. If interpreted as 
an error in the halo mass, a reduction in mass by a factor of 1.6 is required to reconcile the 
model with the upper limi ts, assuming a simple CR luminosity scaling relation, oc M, 



1.46 



(jPinzke &: Pfrommeijl2010l ). or a factor of 4.3 according to Equation lAll in the case when the 
gas density profile is fixed from X-ray observations. For Fornax, the non-detection of a CR 
component is consistent with the low level predicted by the model. 

Assuming no CR emission, we obtain a detection of DM in Virgo at the 4.4(T significance 
level (TS=19.8) for a model in which supersymmetric DM particles of mass w 28 GeV 
annihilate into the 6 h c hannel giving rise to an extended surface brightness profile as predicted 
by the simulations of iGao et al.l (j201ll ) . For Coma and Fornax the significance of a DM 
component is lower, 2.3(7 and 2.1cr respectively. The inclusion of cosmic rays suppresses the 
likelihood of the DM annihilation component, depending on the treatment of the CR emission. 
For Virgo the significance decreases to 3.4(7 when a free CR component is included and to 2.9(7 
and 2.1(7 respectively when our fiducial and optimal CR models are included. The predicted 
DM annihilation flux and cross section upper limits are more robust to the uncertainty in 
the CR component for all three clusters, with the results for different CR models agreeing to 
within a factor of two. 



Assuming neither CR nor DM annihilation emission results in excess emission within three 
degrees around each cluster. The spectrum within this region shows that this excess peaks 
at around 1 GeV. These two features combined lead to an increase in likelihood when an 
extended emission component with an almost scale-free spectrum is included in the model, 
explaining the significance seen in the CR and DM components. 

The DM model that best fits the spectral excess has a particle mass in the range 20 ~ 60 GeV 
for hh final states or 2 — 10 GeV and > 1 TeV for /i^/i^ final states. These mass ranges are 
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consistent with the con clusions of the recent an alysis of a gamma-ray excess in the direction 
of the Galactic Center (jHooper &: Lindenll201lh . The preferred masses are robust to changes 
in the treatment of cosmic rays. Within these mass ranges, DM emission is preferred by the 
data over CR emission in Virgo and Fornax. 

• Models in which the DM annihilation emission has the extended profile predicted by cosmolog- 
ical simulations have a higher significance and higher flux upper limits than models in which 
this emission is assumed to be a point source. Due to the large boost factor of the order of 1000, 
the cross-section upper limits for extended models are at least 100 times lower than those for 
point source models. Our cross-section constraints are much ti ghter than those from an analy- 
sis of clusters using the 11-month data (jAckermann et al.ll20ldl and are also tighter than those 
from a joint analysis of the Milky Way's dwarf galaxies (IGeriiiger-Sameth &: Koushiappas 



2011 



The Fermi-LAT Collaboration: M. Ackermann et al 



2011 



). Our new limits exclude 
the thermal cross-section for < 100 GeV for bb final states and for < 10 GeV for 
final states. The thermal cross-section, however, can still be reconciled with the data 
by lowering the assu med boost factor fro m subhalos and hence the J factor. Since the boost 
factor, b oc M^-?-^^^ (jSoringel et ahlboosl ^. a cut-off mass of 10 '^Mq, rather than our assumed 
1O~^M0, would be sufficient to increase the cross-section limits by a factor of 3. 



In our analysis we have chosen to adopt the values given in the LAT 2-year point source 
catalogue for the intensity and spectral shape of the point sources lying within the virial radius of 
each cluster (where the DM annihilation and CR emission levels are relatively important). This 
allows for possible corrections to the point source parameters in the presence of the CR and DM 
components, while also avoiding the risk of re-fitting sources lying near the boundary of the data 
region with less accuracy. However, we also tried keeping all the point sources fixed or allowing the 
parameters of all the point sources within the data region to vary during the fitting. We find that 
this freedom in the treatment of the point sources has little impact on the DM model fits. Our 
results are insensitive to the uncertainty in cluster halo masses or concentration parameters. A 
change in halo mass by a factor of 4, or a change in concentration parameter by a factor of 2 would 
only change the overall normalization of the inferred DM annihilation emission by a factor of 2 with 
the profile shape remaining unchanged outs ide Q.QlRmr. We al so checked that the different energy 
cuts assumed in our analysis and in that of Imiang et aD (|201lh has no effect on the derived upp er 



limits. We are able to reproduce the upper limits on the annihilation cross-section of iHuang et al 



(|201ll ) for the test case of the Fornax cluster, after correcting for slightly different assumed boost 



factors. 

The Virgo cluster where we have obtained our most significant DM signal has a complex 
structure. Although we have concluded that the central gamma-ray source in M87 does not interfere 
with our detection, the existence of several neighboring structures (M86/M84, M49, M60, MlOO), 
which show signs of interactions, could result in a complex annihilation profile. Given the high 
significance of our result for Virgo, it is worthwhile trying to improve the model for this cluster as 
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a multi-component structure and to use simulations of merging clusters to improve the theoretical 
annihilation profile. 

While we have shown that Virgo has an exceptionally high signal for DM-like emission, the 
signal-to-noise ratio can be enhan ced by stacking several clusters. Such an analysis was recently 



carried out by iHuang et al.l (1201 ll ). but the signal-to-noise was degraded by the assumption of an 
NFW annihilation profile rather than the extended profile seen in the simulations. They considered 
an extended subhalo-dominated annihilation profile but only for individual clusters, not for the 
stack. Their stacked analysis resulted in looser constraints on DM annihilation than their analysis 
of individual clusters, presumably because the use of an inappropriate theoretical profile resulted 
in the different clusters yielding inconsistent results. Thus, it is clearly worth repeating the joint 
analysis with the "correct" subhalo-dominated profile. It is also tempting to extend the search for 
DM annihilation to multi- wavelength data from the radio to very high energy gamma-rays. 

The significance of the excess DM-type annihilation emission that we have detected in our small 
cluster sample is only marginal and needs to be confirmed with future observations. Nevertheless, 
the similarity with the excess gamma-ray emission from the Galactic center and the general trend 
for such excess emission in all three clusters examined here is intriguing. 
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A. Semi- Analytic formula for the Cosmic Ray induced gamma-ray emission 



Here we summarize the relev ant equations for calc u lating the CR induced ganima-r ay emission 
in galaxy clusters as derived in iPinzke &: Pfrommerl (j2010l ) and iPinzke et al.l (j201ll ). The CR 
induced photon source function from pion decay can be decomposed as: 



dtdVdE 



A{r)s{E). 



The spatial part is given by: 

A{r) = aa 



vir Ccenter){^ ~^ ^ D ) ^) ~^ Ccenter) Pgasi'^} 

^trans 



with 



Rtrans 



The spectrum is given as: 



1.7 X lO-'^ X (M^ir / 10^^ Mq)^-^^ 
-- OmiR^ir X (M^,J1O15M0)O-39 

= 1.04 X {A'Uir/W^^MQf-^^ 



t ■jrVJ _j-U 



with A = (0.767,0.143,0.0975), a = (2.55,2.3,2.15), nip the proton mass and 

E^ 



Dy{Ey, Ey^break) — [1 + ( 



^7 ^j3i-l/9 



The proton cut-off energy is 



10« 



-^■y, break 

R 



and 



app^i ~ 32(0.96 + e^-''2-2-''"0mbarn. 



The maximum shock acceleration efficiency is c hosen to be Cn 



gas density is fitted with multiple beta-profiles (jPinzke et al 



max. 



2011 



0.5 so that g{Cp,r 



(Al) 



(A2) 
(A3) 
(A4) 



as: 



Pgas 



nir, 



XhX( 



1/2 



(A5) 

(A6) 

(A7) 

(A8) 
1. The 

(A9) 



where Xh = 0.76 is the primordial hydrogen mass fraction and X^ = 1.157 is the ratio of electron 
and hydrogen number densities in fully io nized ICM, with parameter values for nj(0), r^^i and /3i 
listed in TABLE VI of lPinzke et"aD (|201lh . 



